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Abstract: We present a general methodology for performing statistical infer- 
ence on the components of a real- valued matrix parameter for which rows and 
columns are subject to order restrictions. The proposed estimation procedure 
is based on an iterative algorithm developed by Dykstra and Robertson (1982) 
for simple order restriction on rows and columns of a matrix. For any order re- 
strictions on rows and columns of a matrix, sufficient conditions are derived for 
the algorithm to converge in a single application of row and column operations. 
The new algorithm is applicable to a broad collection of order restrictions. In 
practice, it is easy to design a study such that the sufficient conditions derived 
in this paper are satisfied. For instance, the sufficient conditions are satisfied 
in a balanced design. Using the estimation procedure developed in this article, 
a bootstrap test for order restrictions on rows and columns of a matrix is pro- 
posed. Computer simulations for ordinal data were performed to compare the 
proposed test with some existing test procedures in terms of size and power. 
The new methodology is illustrated by applying it to a set of ordinal data 
obtained from a toxicological study. 

1. Introduction 

In many applications the parameter of interest 9 can be expressed as elements of a 
real- valued I x J matrix such that the elements within rows and/or columns of the 
matrix are subject to inequality restrictions called order restrictions. Researchers 
are often interested in drawing statistical inferences on such parameters subject to 
a variety of order restrictions. For a parameter vector rj = (?yi , 772 , . . . ,r] p )', some 
commonly encountered order restrictions are; simple order, where rji < 772 < • • • < 
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Table 1 

For each response variable, the data structure of the experiment in Wormser et al. [26] 



Genotype 


Sample 
size 




Level of skin injury 




Unremarkable 


Minimal 


Mild 


Moderate 


Marked 


Cox-2 deficient 


10 


"■1,1 


"1,2 


"1,3 


"1,4 


"1,5 


Wild type 


10 


"2,1 


"2,2 


"2,3 


"2,4 


"2,5 


Cox-1 deficient 


10 


"3,1 


"3,2 


"3,3 


"3,4 


"3,5 



•q p , umbrella order, where 771 < 772 < • • • rji > iji+i > • • • > rj p and simple tree order ; 
where 771 < r]i, for all 2 < i < p. 

Recently Wormser et al. [26] conducted an experiment to evaluate the differ- 
ences among three different genotypes of mice, namely, the wild type (WT), the 
cyclooxygenase-1 deficient (COX-l-d) and the cyclooxygenase-2 deficient (COX-2- 
d) mice, when they were exposed to sulfur mustard (also known as mustard gas). 
Depending upon the severity of injury to skin, each animal was categorized into 
one of five ordered categories, namely, "unremarkable", "minimal", "mild", "mod- 
erate", and "marked" (details in Section 5). In this experiment, a sample of n = 10 
animals from each genotype was exposed to sulfur mustard. Let n^j, i = 1,2,3, 
j = 1, 2, . . . , 5, denote the number of animals in the i th genotype that belong to 
jth reS p 0nse category, with E(riij) = mrij. Then the parameters of interest are 
the cumulative probabilities 9i_j = J2k=i n i,k> * = 1,2,3, j = 1,2,3,4. Note that 
6i t s = 1, i = 1,2,3. Table 1 summarizes the type of data obtained in Wormser et 
al. [26]. Clearly, the rows of 9 satisfy a simple order as they represent cumulative 
probabilities. According to Wormser et al. [26], COX-2 deficiency has a protective 
effect against inflammatory processes while COX-1 deficiency has a negative effect. 
Consequently, each column of 9 is also subject to simple order restriction. Thus in 
this example the rows as well as columns of 9 are subject to simple order restriction. 

The above type of matrix order restrictions commonly arise in a variety contexts 
such as the analysis of ordinal data (cf. Agresti and Coull [1], Grove [9], Nair [16] 
and Wang [25]), survival analysis (Praestgaard [19]), Phase I clinical trials involving 
"cocktail" of treatments (Conaway et al. [6]) and analysis of time-course and dose- 
response gene expression microarray studies, etc. 

For a given parametric family with matrix valued parameter 9 G O C 1Z IxJ , 
where O is the parameter space defined by order restrictions on the rows and/or 
columns of 9, one may estimate 9 using restricted maximum likelihood estimators 
(RMLEs) and test alternative hypotheses using the likelihood ratio tests and their 
modifications. Such methods have been well studied in the literature (cf. Barlow 
et al. [2], Robertson et al. [21] and Silvapulle and Sen [23]). However, as seen from 
Hwang and Peddada [11] and Lee [12], the RMLE is not always efficient relative to 
the unrestricted maximum likelihood estimator (UMLE). Also, the likelihood ratio 
tests in the present context may not be computationally simple and the asymptotic 
distribution under the null hypothesis may involve nuisance parameters (cf. Franck 
[S], Grove [9], Robertson and Wright [20] and Wang [25]). Silvapulle [22] provided 
an interesting explanation for why sometimes the RMLEs and the likelihood ratio 
tests provide counter-intuitive results. 

In view of the general concerns regarding RMLE and the likelihood ratio tests, in 
Section 2 we introduce a computationally straightforward methodology for estimat- 
ing a matrix valued parameter 9 when the rows and/or columns are subject to order 
restrictions. The proposed estimation procedure is based on the iterative procedure 
of Dykstra and Robertson [7] and uses the method of Hwang and Peddada [1 1] for 
general order restrictions. If the elements of the random matrix are independently 



61 



E. Teoh et al. 



and normally distributed and if rows as well as columns of the matrix are subject 
to simple order restriction then the proposed procedure is identical to the estimator 
given in Dykstra and Robertson [7], but the two procedures may differ for other 
order restrictions. Using the proposed point estimators, in Section 3 we introduce 
a Kolmogorov-Smirnov type test statistic and a bootstrap-based methodology for 
determining significance. The performance of the proposed test procedure is eval- 
uated using computer simulations in Section 4 and in Section 5 it is illustrated by 
analyzing the data in Wormser et al. [2G] mentioned above. Concluding remarks 
are provided in Section 6. 

2. Estimation of parameters subject to order restrictions 
2. 1 . Notations and a brief review 

Throughout this paper 1Z P denotes the vector space of p x 1 real vectors and TZ IxJ 
denotes the vector space of / x J real matrices. Two components 9ij and 9 T:S of 
9 £ C 1Z IxJ are said to be linked if the inequality between them is known a priori. 
A parameter is said to be a nodal parameter if it is linked to all IJ components of 
9 £ O C 7Z IxJ . A subset of parameters A4, formed by the components of 9 S C 
1Z IxJ , is a linked subgraph if all parameters in M. are linked, with at least one strict 
inequality. Note that every linked subgraph represents a simple order-restriction 
and conversely, every simple order is a linked subgraph. A linked subgraph M. is 
said to be maximally linked if for any linked subgraph M , A4 C J\f ==>■ M. ~ M . 

If a subset of parameters simultaneously satisfies two linked subgraphs M. and 
M , then we use the notation MAN to describe the subset. As noted in Peddada 
et al. [18], any order-restriction between parameters can be expressed in terms of 
a collection of maximally linked subgraphs. Similarly, every parameter 9ij appears 
in a finite collection of maximally linked subgraphs as a nodal parameter (within 
each subgraph). Such maximally linked subgraphs are said to be associated with 

If the rows of 9 are subject to order restriction 3? C 1Z 1 and the columns are 
subject to order restriction C C 1Z J then we shall denote the joint order restriction 
by RAG C K lxj . 

The notation 31A1Z J would indicate order restrictions 31 on the rows and no order 
restrictions on the columns of 9. Similarly, the notation 1Z 1 AC would indicate order 
restrictions C on the columns and no order restrictions on the rows of 9. 

It is important to emphasize that in this article we only consider linked subgraphs 
which are subsets of rows or which are subsets of columns of the parameter 9. Thus 
we are not considering inequalities between arbitrary linear or non-linear functions 
of the elements of 9. 

Example 2.1 (Umbrella order restriction on the rows of 9). Suppose the ele- 
ments of each row of 9 are subject to an umbrella order with peak in the s th 
column. That is, the components of the i th row of 9 satisfy the inequalities 6^1 < 
9i.2 < •• < 9i. s > #j jS +i > ••• > Oi ; j, i = 1,2,..., I. Then in this case the 
parameters of the i th row can be expressed in terms of two maximally linked 
subgraphs, namely, M^-.s = #i,2, • • • I < #i,2 < ••• < #m} and 

Mi^-.j = {(9i, s ,6 iiS+ x,...,6 it j) | 6> M > 6> M+ i > ■■• > Qi,j}. Thus, the subset of 
parameters in the i th row of 9, i = 1,2, ... ,1, can be expressed as jVfi,i :s AMi. s:i /. 
Further, for r < s, the maximally linked subgraph associated with 9i :T in the i th 
row, i — 1,2,...,/, is M.i,x-. s . 



Order restrictions on rows and columns of a matrix 



65 



For a vector x = (x\, x<i, . . . , x p )' £ 1Z P , a simple order operator £® : 1Z P — ► S is 
an orthogonal projection operator onto the simple order cone § = {i£ 1Z P : x\ < 
x 2 < ■ ■ ■ < x p } where w is a vector of positive weights and the i th component of 
is given by 

(2.1) = minmax ^' 3 3 



s>i t<i T,j=t w j 

We shall use the terminology simple order function to describe the min-max formula 
used in (2.1). 

Remark 2.1. For a fixed weight vector w, is a monotonic operator in the sense 
that if x < y then £%{x) < €f u (y), where the inequality is componentwise. 

When estimating a parameter rj = (r/i, 772, . . . , r} p )' subject to arbitrary order 
restrictions 6 C W with at least one nodal parameter, Hwang and Peddada [11] 
used the simple order operator with a suitable weight vector w, for estimating 
the nodal parameters. Typically the weights are proportional to the precision (or 
sometimes the sample size) of UMLE. Once a parameter is estimated, then the 
corresponding UMLE is replaced by the new restricted estimator and is assigned 
arbitrarily large weight B, B — > 00, in all subsequent calculations. To estimate 
a non-nodal parameter, identify the collection of all maximally linked subgraphs 
associated with that non-nodal parameter and apply the simple order operator 
on the vector corresponding to the subgraph with a suitable weight vector w. 
Suitable modifications were proposed for graphs with no nodal parameters. 

Example 2.2 (Umbrella order restriction). Suppose rj is a parameter satisfying 
the order restriction 771 < 772 < T]3 > T}a > V5 with UMLE fj. In this case the 
only nodal parameter is 773. Suppose w = (wx, W2, W3, W4,, W5)' is the weight vector 
associated with fj. According to Hwang and Peddada [11] the estimation procedure 
begins with the nodal parameter 773. For i = 1,2, denote wn\ = Wi, f)a-\ = i)i and let 
iU( 3 ) = w 5 , i/J( 4 ) = 7774, W(5) = w 3 and r} (3) = 775, fj^ = 774, fy 5) = i) 3 . Then 773 may be 
estimated by the following simple order formula 

773 = max— 7 . 

Next, the non- nodal parameters 771, 772, 774 and 775 are estimated using the maximally 
linked subgraphs 771 < 772 < 773 and 773 > 774 > 775, respectively. The estimators for 
the non- nodal parameters 771,772 are simplified as follows: 

i . r „ Wifji + 77J 2 772 x , ; . . ,W\f\\ + w 2 rj 2 . , s , 

771 = mm{77i, ■ ,773}, 772 = mm{max( ■ , ?72 ) , ?73 1 - 

Wi +w 2 Wi + W2 

In a similar manner 774 and 775 are estimated. 

Remark 2.2. For a given w, €^ is a function of several simple order functions and 
therefore £^ is a monotonic operator. That is, for all x < y, ^(x) < <£^(y), where 
the inequalities are componentwise. 

2. 2. Estimation of matrix valued parameters 

We extend the notations from the previous section to matrix valued parameters 
as follows. For a matrix X £ 1Z IxJ and a weight matrix W, we denote the matrix 
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simple order column operatorhy C^. Each column x of X is orthogonally projected 
onto the simple order cone § using (£® , where w is a suitable column vector of W . 
Similarly, the matrix simple order row operator \s denoted by ytfy, which, with rows 
of W, projects each row vector of matrix X £ 1Z IxJ orthogonally onto the simple 
order cone S. Analogously, for arbitrary order restrictions CD, we define matrix 
column and row operators using Hwang and Pcddada methodology by and 
91^, respectively. 

We now describe the proposed algorithm for estimating 9 £ = CRAC C 1Z IxJ 
using an unrestricted point estimator 9. We use the weight matrix Wr when op- 
crating on the rows of 9 and weight matrix Wc when operating on the columns 
of 9. 

Step 1 (An unrestricted estimator): 

Obtain an unrestricted estimator 9 for an / x J matrix parameter 9. Note that in 
most situations a user may prefer to start with the unrestricted maximum likelihood 
estimator (UMLE), although it is not required. 

Step 2 (Estimation under order restrictions on the columns of 9): 
Apply the procedure of Hwang and Peddada [] 1] on each column of 9 to obtain 
estimates for 9 under the order-restriction on the columns of 9. That is, apply the 
operator on the columns of 9 and denote the resulting estimator by €.^ Vc (9). 
Note that the elements of (ff) may not satisfy the order-restriction on the rows 
of 6. 

Step 3 (Estimation under order restrictions on the rows of 9): 

Apply the operator SH^-„ on the rows of (8) and denote the resulting estimator 

by(*^ o <L% o ){6). R 

Step 4 (Iterate to convergence): 

Repeat Steps 2 and 3 to obtain the q th iterate, ffiw R ° £h/ c )'(^)- Stop when some 
reasonable convergence criterion is reached. Denote lim^oo (^w R ^w^ c ) 9 (^) as 
k- 

Step 5 (Final estimate): 

Since 9t^- R and ^%/ c do not necessarily commute for all order restrictions and 
weights used in the calculated weighted averages, repeat the process beginning with 
within-row order restrictions followed by within-column order restrictions. That is, 
compute 9i = lim^oo (<£^ c ° ^w R ) q (^)- The final estimate of 9 £ is taken to 
be 

(2.2) J = I(0 1 + 2 ). 

Remark 2.3. In general 9\ 7^ 02- To illustrate this, consider the very special case 
where 9 is a 2 x 2 matrix with rows and columns both subject to the simple order 
restriction 9ij < 6*2j, < 6^2, i = 1,2, and j = 1,2. Consider the extreme case 
where 9 is a symmetric matrix with 9\ t \ < #1.2 = 82,1 > $2,2- Further, suppose 
Wc = Wr = 11'. Thus we have a perfectly "symmetric" problem where we may 
expect 6*i = 9i- However, even in this rather seemingly obvious situation 9\ 7^ 62, 
but 9\ = 9' 2 . Hence the composition ° ff|r c ) is not commutative. For this 

reason, we need to invoke Step 5 in all situations. 

Before we discuss the convergence of the above algorithm in Steps 4 and 5, we 
consider the following example which may be instructive. 
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Example 2.3. Consider a clinical trial comparing 3 new treatments with an ex- 
isting treatment using 4 dose groups per treatment group. For the j th dose of the 
i th treatment, i, j = 1, 2, ... ,4, let 9ij denote the sample mean response based on 
riij = n observations. Further, for simplicity of illustration, let 9ij ~ mde P N(9i i j, c), 
i, j = 1, 2, 3, 4. Without loss of generality, let the first row of 9 correspond to the ex- 
isting treatment. Then the order restrictions of interest are 6\j < 6ij, i > 2, j > 1, 
i.e. simple tree order within each column of 9, and 9ij 1 < 9\ J2 , 1 < ji < j'2 _ 4, 
i > 1, i.e. simple order within each row of 9. We choose Wc = Wr = 11', where 
1 = (1,1, 1,1)'. 

We begin with Step 2 of the algorithm. Let Y = €.^ Vc {9). Thus columns of Y 
satisfy simple tree order. That is, 

Yij>Y ld ,Vi = 2,3,4, j = l,2,3,4. 

Now applying Step 3 on Y = &%r {9) wc obtain Z = ^ R (Y). The rows of Z 
satisfy a simple order. That is, 

Zi,j t > Z ith , V ji,j 2 = 1, 2, 3, 4, ji > j 2 , i = 1, 2, 3, 4. 

We now demonstrate that the columns of Z would also satisfy the simple tree order 
restriction imposed on 9. That is, for any column j we need to demonstrate that 
Zi,j > Z\ j , for all i = 2, 3, 4. Note that for each j = 1, 2, 3, 4, 

, = minmax ^^ ,Vi = 1,2,3,4. 
J t>j «<i t - a + 1 ' 

Since F.fc > F.fejVi = 2,3,4, = 1,2,3,4 and since the above simple order 
function is a monotonic function it follows that for all i = 2, 3, 4 and j = 1, 2, 3, 4, 

Z t j = mm max — — > mm max — = Z\ ,• . 

*>i s<j t - s + 1 t>j s <j t - s + 1 ,J 



o 



Thus 9i = (Pvw R o £{y )(£)• Similarly, it can be demonstrated that 9i — (^w c 

)W' Thus in this example the algorithm converges after one application of 
column and row operations and q = 1. 

As will be demonstrated formally in the following theorem, one of the reasons 
for the convergence observed in the above example is that Wc = Wr = 11', where 
1 is a column vector of l's of suitable length. Or more generally, Wc and Wr 
are each of rank 1. In many applications, researchers use a balanced design for 
all dose and treatment combinations. In such situations it is appropriate to take 

W C = Wr = 11'. 

We now discuss convergence of Steps 4 and 5 in the above algorithm in the 
following theorem. 

Theorem 2.1. Suppose 9 6 = C 1Z IxJ , with every row subject to the same 
order restriction and every column subject to the same order restriction. However, 
the order restriction on a row need not be same as the order restriction on a column. 
Further, suppose that the weight matrices Wr and Wc are each of rank 1. Then 
for all X G TZ IxJ , 

(a) € e Wc o (X* R o tfr a )(X) = (X* H o € e Wc )(X) e 9. 

(6) ml R o (<^ c o m* R )(x) = {& Wc o mlj(x) e e. 
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Thus it takes one column operation and one row operation for the algorithm de- 
scribed in Steps 4 o- n d 5 to converge. 

Proof. We prove the theorem for (a) since the proof of (b) follows similarly. The 
main underlying idea of the proof is that, as stated in Remark 2.2, the row and 
column operators an d ^w c are monotonic and the simple order function used 
in these operators is a function of suitable weighted averages of the elements of #. 

Let Y = €y t , c {X). Thus, the columns of Y satisfy the order restriction 6. That 
is, for each j = 1, 2, . . . , J, and any two linked parameters 9 ix j and #j 2 j, with 
lu3 < e i2tj ,we have Y iuj < Y 42J , Let Z = tf* R (Y) = o' <L% C ){X). Thus, 

the rows of Z satisfy the order restriction "R. 

For each j = 1,2,..., J, for any two linked parameters and #i 2j , with 
9i 1 ,j < 9i 2 j, we demonstrate that Zi lt j < Zi 2 j. 

Recall from Remark 2.2 that the operators and 9^w R are functions of several 
simple order functions of suitable weighted averages of the components of 6. For a 
given row i\, let C denote the set of all subsets of J = {1, 2, 3, . . . , J} such that 
8 with column indices in these sets are used in the construction of j. Further, 
since the order restrictions in every column is the same, the same set of column 
indices are used in the construction of Zi 2 j. Thus Zi u j and Zi 2 ,j are functions of 
J2keK( W R)n,kY luk /J2keK( W R)n,k and Y, keK {W R )i 2 , k Y i2>k / T, keK (W R ) i2 , k , re- 
spectively, where K C C. 

Since 9^w R 1S a monotonic operator, it is sufficient to prove that 

( 3) Z^WnKk ~ E k&K (w R ) i2 , k ' e c 

Since Wr is of rank 1, we can therefore express (Wr)^^ — ca 2 {Wii)i 2 ,k- There- 
fore 



(2.4) 



J2keK( W R)n,kYt uk _ Y^keK a i^R)i2,k Y ii,k _ Efeeif(W /r fl)i2,fc^i,fe 



But since j < j for every j — 1, 2, . . . , «/, (2.4) is bounded by 

T,keK( W R)i2,kY i2:k 
J2keK( w n)i2,k 

Thus, by the monotonicity of the operator Zix,j < Zi 2 j. Hence ffiw ° 

£^)(I) € 9 and therefore <L% c is a left identity of {^ Vr o €^ c )(X). 

Hence the proof of the theorem. □ 

Remark 2.4. An important consequence of the above theorem is that Wc and Wr 
need not be limited to the matrix 11', but Ix J weight matrices of the type nij = nj, 
i = 1,2,...,/, j = 1, 2, . . . , J, can be considered. For example, in a treatment by 
dose-response study, it is not required to have equal sample sizes in all treatment 
by dose combinations, but the experimenter may conduct a study with a sample 
of n.i for the i th treatment, as long as within each treatment the same number of 
observations are collected on each dose group or vise versa. 

Remark 2.5. As can be seen from the following counter example, it is encouraging 
to note that the sufficient condition stated in the above theorem is not necessary. 
Let 9 be a 2 x 2 matrix with 6^1 < #i,2, i = 1,2, and 9\,j < #2j, j = T 2. Suppose 
that #ii < #i2, #i.i < #2.1 < #2 2, an d #1,2 > #2,2- I n this case, for any non-singular 
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weight matrices Wr and Wc, the iterative algorithm converges in one cycle, i.e. 
one column operation followed by one row operation. 

Remark 2.6. In some situations, such as in the context of ordinal data, the rows 
(or columns) of the unrestricted estimator 9 naturally satisfy the order restriction. 
In such cases, under the conditions of Theorem 2.1, we need to apply on 9 (or 
5Rff R ) only once. 



2.3. A simulation study 

For an Ix J matrix 9 whose components are independently normally distributed, we 

compared the performance of the proposed order-restricted matrix estimator 9 with 
the UMLE 9 using a simulation study. Although we considered a variety of order 
restrictions, patterns of means, sample sizes and weight matrices, in this article we 
provide the results of a small subset since very similar results were obtained across 
all patterns. 

In the simulation study reported here, E(9) = 9 with 9i.j = i + j, 1 < i < I, 
1 < j < J - 2, and 9 l . J = i - j + J + 1, 1 < i < I, j > J - 2. Thus we have 
a simple order along the columns of 9, with 9i x j < 9i 2 j, for all 1 < i\ < %2 < I, 
j = l,2,..., J, and a tree order along the rows of 9, with Oi l < 6ij, for all 1 < j < J 
and 1 < i < I. Each of the normal variables was generated with a standard deviation 
of 1 and we chose sample sizes Rj j = i, i = 1, 2, . . . , I, and j = 1, 2, . . . , J so that 
Variance (6i,j) = l/i- The (i,j) th element of the weight matrix Wr is given by \fi 
and the weight matrix Wc = W' R . 

Our simulation study is based on 10,000 simulation runs and the results are 

summarized in Table 2. In addition to comparing the average bias of 9 with that of 

9, we also computed the percentage reduction in quadratic and quartic loss due to 

s „ V'_ V J _ E&j-eij)* 
9. The reduction in loss relative to 9 is defined as 100 x (1— _ j 1 ,_ J ; 1 — ; ), 

where 6 = 2 corresponds to quadratic loss and 6 = 4 corresponds to quartic loss. 

Observe that the proposed procedure reduces the average quadratic loss and 
quartic loss substantially, without costing much in terms of bias. 



3. Testing hypotheses under order restrictions 

In some applications, such as in Wormser et al. [26], researchers are interested in 
performing tests of hypotheses regarding the elements of each column of 9 when 
the rows are subject to order restrictions. The order restrictions on the rows are 
not part of the hypothesis, but they are present due to the underlying probability 
model or for other reasons. Thus, the hypotheses of interest are 

Hq '■ &x,j = 92 j = ■ ■ ■ = Ojj, 1 < j < J, 



Table 2 



Bias and reduction in loss due to the proposed estimator 8 relative to UMLE 8 



I 


J 


Average bias 


Percentage reduction 


relative to 












Quadratic loss 


Quartic loss 


2 


5 


-0.0774 


0.0014 


29.80 


53.54 


2 


10 


-0.0608 


-0.0009 


22.38 


41.41 


5 


5 


-0.0612 


0.0018 


28.12 


55.06 


5 


10 


-0.0324 


-0.0011 


22.75 


45.02 
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(3.1) H a : (e hj , e 2 j, . . . , e Itj y € A p k=1 M k , 1 < j < J, 

where Aik are maximally linked subgraphs. However, each row of 9 is itself subject 
to the restriction (6i i, 0i.2, ■ ■ ■ , Qi.j)' € A q k=1 Mk, ^ — * — where are maximally 
linked subgraphs. 

In other examples researchers may be interested in testing 

Ho : 9 itj = B V j. V(i, j) + (if,?), 1 < < 1, 1 < j, j' < J, 



(3.2) iJ a : G IRAC, 

where C = Aj^A'ffc and 3? = A^ =1 A/fc, and each M.u and A/fc are suitable maximally 
linked subgraphs. 

In both (3.1) and (3.2) the point estimators of 9 are the same, and are obtained 
under the order restrictions KAC using the methodology introduced in Section 2, 
although the test statistics are different. 

For a maximally linked subgraph Mk, defined by 9 s _k < S s +i,fe • • • < 9 r _k, the 
two parameters 9 s- k and 9 r _k, which arc at the ends of the graph, are said to be 
farthest linked parameters of the subgraph. Under this maximally linked subgraph 

let 9 St k, s +i,fc ■ ■ ■ , 9 r .k denote the estimated value of the parameter (9 S} k, #s+i,fc, ■ ■ ■ , 
9 r ,k)' using the methodology described in Section 2. To test the hypotheses given 
in (3.1), for each maximally linked subgraph, compute the estimated difference 
between the two farthest linked parameters of the subgraph and divide it by the 
standard error of the difference under the null hypothesis and take the largest 
over all maximally linked subgraphs. More precisely we propose the following test 
statistic for testing (3.1): 

a , _ a , 

(3.3) Xi = max 



where 9 r ^ and 9 S ^ are the farthest linked parameters in the maximally linked 
subgraph M.^. and the max is taken over all maximally linked subgraphs Mk- 
Similarly, we may test 

Hq : 9i t i = 9i t 2 = ■ ■ ■ = 9i j, 1 < i < I, 



(3.4) H a : 6 ii2 , ■ • ■ , 9 hJ y G A£ =1 A/" fe , 1 < i < I, 
using the test statistic 

(3.5) T 2 = max 6k f ~ 0k f 

Aft se{6 k>r - 6 k , s ) 

where 9k, r and 9k. s are the farthest linked parameters in the maximally linked 
subgraph A4 and the max is taken over all maximally linked subgraphs A/fc . Then 
the hypothesis on rows and columns (3.2) may be tested using T = max{Xi, T2}. In 
the above expressions se is computed under H . A suitable such estimate is derived 
in Section 4 for the case of ordinal data. The critical values and p- values are obtained 
by the bootstrap methodology as follows. For each bootstrapped dataset, the i th 
independent group is formed by taking a simple random sample (with replacement), 
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of appropriate size, from the pooled sample of all subjects across all independent 
groups. In this resampling procedure, we sample the entire record of a given subject. 
For each bootstrap sample we compute the test statistic T, which is denoted by 
T* . The sampling distribution of T* is obtained by generating a large number 
of bootstrap samples, say 10,000. Then the bootstrap p-value is computed as the 
proportion of times T* exceeded the observed T. 

4. Simulation study 

We conducted a simulation study to investigate the performance of the proposed 
bootstrap test based on T\ for comparing I treatment groups when the responses 
are measured on J ordered categories. In a random sample of n observations on 
the i th treatment, i = 1,2,..., I, let X^j denote the frequency of responses in 
the j th ordered category, j — 1,2,..., J. Further, let E(Xij) = nxjj and let 
&i>j = Si=i 7r *-fc denote the cumulative probabilities with O^j = 1. Let -fr^.j denote 
the sample proportion Xij/n and let 6ij denote the corresponding cumulative sum 
of sample proportions. Thus under the multinomial model 9 is the UMLE of 9. Note 
that the rows of 7r, and hence the rows of 9, are independently distributed. 

We performed simulations to study the size and power of T\ when testing Hq 
against H a — H , where H : 9 e O = {9 \ 6 r j = 9 s j, 1 < r, s < I, j = 
1,2,..., J} and H a : 9eQ = {9\r<s^> 9 r j < s j}. Ordinal data was 
generated for a variety of parameter configurations (Tabic 3) with sample sizes of 
10, 20, and 50 subjects for each independent group. 

Under the null hypothesis, we estimate the standard error of (9 r _j — 9 s _j) re- 
quired in Ti by s"e(9 r j - § Stj ) = \f^j, where V 3 = \ J2i=i Z)»=i[*r(l ~ K r )I(r = 

s) — TT r -if s I(r 7^ s)], and 7T r = ~ ; r = 1, 2, . . . , J, a pooled Bayes esti- 

mator under quadratic loss and suitable Dirichlct distribution prior (Lehmann [15], 
page 293). Since the usual pooled MLE can take a value of or 1 with a positive 
probability, we prefer the above estimator over the MLE for calculating Vj. 

For each configuration listed in Table 3, we performed 10,000 simulations to 
evaluate the size and power of T\. Critical values of T\ were determined using 
10,000 bootstrap samples for each simulation run. The nominal size was set at 
a = 0.05. We compared the proposed test with the order-restricted methods of 
Grove [9] and Nair [16] as well as with the one-sided Kolmogorov-Smirnov test. For 
simplicity, we bootstrapped the critical values of Kolmogorov-Smirnov test. Since 
the procedure of Grove [9] is designed for comparing only two groups, we compared 
the performance our procedure with Grove [9] for / = 2 only. Similarly, comparisons 
with the one-sided Kolmogorov-Smirnov test was limited to the case 7 = 2 only. 

Results of the simulation study are summarized graphically in Figure 1 using 
scatter plots. The top three panels correspond to the type 1 errors for the three 
different sample sizes per test group, i.e., n = 10,20,50. The bottom three panels 
correspond to the power for the corresponding sample sizes. In each panel the ver- 
tical axis represents the estimated probability of rejection of null hypothesis by the 
proposed critical region, while the horizontal axis represents the estimated proba- 
bility of rejection of null hypothesis by the three alternative procedures. Thus the 
six scatter plots represent the comparison between the proposed and each of the 
three competing test procedures. The '+' symbol corresponds to the comparison 
between the proposed test and Grove's test (Grove [9]), corresponds to compar- 
ison between the proposed test and Nair's test (Nair [16]) and '*' corresponds to 
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Table 3 

Multinomial parameter configurations used in the simulation study 



Group 1 Group 2 Group 3 

Til 7r 12 ^13 7T 14 7T 15 7T 2 l 7T 2 2 ^23 7r 24 ^25 "31 7r 32 7r 34 135 



H 1 = 2, J 



2, J 



: 3, J 



7" = 3, J 
Hi -Ho 1 = 2, J 

1 = 2, J - 



3, J 



1=3, J=5 
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the comparison between the proposed test and the one-sided Kolmogorov-Smirnov 
test. In each scatter plot the diagonal line corresponds to the line of equality. Ad- 
ditionally, in the top three panels a horizontal and a vertical line is provided at 
0.05 + a/(.05 x .95)/10000 to indicate points that exceed the nominal value of 0.05 
by at least one standard error. 

We notice that in general all three test procedures approximately maintain the 
nominal size of 0.05. In situations involving rare events (e.g. - probability vector 
(0.01, 0.01, 0.98) all tests are conservative, but even in that case, relative to others, 
the proposed test appears to recover the nominal size more quickly as the sample 
size increases. 

As indicated by the points above the diagonal line in the three bottom panels 
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Fig 1. Power and size comparisons of the proposed procedure with Grove's test (+), N air's test 
(#) and Kolmogorov-Smirnov test (*). Nominal size is 0.05. Results for the proposed method are 
plotted on the vertical axis, and results for other methods are plotted on the horizontal axis. 



of Figure 1, the proposed test seems to enjoy higher power than its competitors 
in almost all situations considered in this simulation study. In particular, even for 
parameter configurations that are very close to the null, the proposed procedure 
has higher power than its competitors. Additionally, in these cases, the proposed 
procedure appears to increase in power with sample size at a faster rate than its 
competitors. 

In addition to a gain in power, a distinct advantage of the proposed test over 
likelihood ratio type procedures is the ease of implementation for any arbitrary 
order restriction on the rows and columns. The procedure described in Grove [9] is 
limited to two groups. As seen in Wang [25], the likelihood ratio type procedures 
are very challenging to implement as the number of groups increases. 

5. Illustration 

In the experiment of Wormser et al. [26] mentioned in the introduction of this 
paper, the effect of mustard gas on the skin of mice was evaluated using 6 ordi- 
nal variables, namely, subepidermal microblister formation, epidermal ulceration, 
epidermal necrosis, acute inflammation, hemorrhage, and dermal necrosis. The ex- 
periment consisted of exposing 10 mice of each genotype (i.e. COX-2-d, WT and 
COX-l-d) to 0.317mg of sulfur mustard. Changes in skin condition of each animal 
(as measured by the above 6 variables) were noted on ordinal scale ranging from 
"unremarkable" , "minimal" , "mild" , "moderate" , to "marked" . For each response 
variable and each genotype, in Table 4 we provide the sample cumulative proportion 
of animals in each category. 
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Table 4 

Cumulative relative frequencies for each genotype each for response variable 



Level of skin injury 



Response variable 


Genotype 


Unremarkable 


Minimal 


Mild 


Moderate 


Marked 


Microblister 


COX-l-d 





0.2 


1 


1 


1 




WT 





0.5 


0.9 


1 


1 




COX-2-d 


0.3 


O.S 


1 


1 


1 


Ulceration 


COX-l-d 


0. 1 


0.5 


0.8 


0.9 


1 




WT 


0.8 


0.9 


1 


1 


1 




COX-2-d 


1 


1 


1 


1 


1 


Epidermal necrosis 


COX-l-d 











0.3 


1 




WT 








0.1 


0.8 


1 




COX-2-d 


0.1 


0.2 


0.7 


0.8 


1 


Acute inflammation 


COX-l-d 








0.5 


1 


1 




WT 








0.6 


1 


1 




COX-2-d 


0.1 


0.5 


1 


1 


1 


Hemorrhage 


COX-l-d 





0.1 


0.9 


1 


1 




WT 








1 


1 


1 




COX-2-d 


0.1 


0.1 


1 


1 


1 


Dermal necrosis 


COX-l-d 





0.1 


0.8 


1 


1 




WT 





0.1 


0.9 


1 


1 




COX-2-d 


0.2 


0.4 


1 


1 


1 



For each response variable, the statistical hypothesis of interest was motivated by 
the underlying biology. COX-2 is involved in a variety of inflammatory processes 
caused by noxious stimuli. For instance, COX-2 is induced within 12 hours and 
persists up to 3 days after excisional injury in rat skin. COX-2 expression is seen 
in the basal cell layer, peripheral cells in the outer root sheath of hair follicles, and 
in fibroblast-like cells and capillaries near the wound edges Lee et al. [13]. 

Neutrophil COX-2 protein expression after burn-induced injury in mice signifi- 
cantly increased at 4 hours and dramatically decreased 36 hours after injury (He 
et al. [10]). Due to the central role of COX-2 in the inflammatory processes, it is 
believed that the effect of mustard gas on COX-2 deficient animals will tend to be 
less severe than that on the Wild Type animals. 

Nevertheless, COX-1 may have a protective effect on the skin against sulfur 
mustard as shown in other organs like the kidney and brain (Lin et al. [14]) and 
Vane et al. [24]. COX-1 was suggested to confer protection on the epithelial cells 
of the crypts of Lieberkuhn, through promotion of crypt stem cell survival and 
proliferation, in the ileum of irradiated mice Cohn et al. [5] . 

As in our skin model case, Cohn et al. [5] also concluded that prostaglandins pro- 
duced through the COX-1 pathway, may not be important in unstressed conditions, 
but still may have a protective role in the response to epithelial injury. Therefore, 
we expect mice with COX-1 deficiency to have a more severe response, on average, 
than wild- type animals. Then, the experimental setting is the comparison of three 
groups with six responses, each of which is measured on an ordinal scale with five 
levels and subject to a simple order. 

We analyzed each response variable separately but adjusted the p-value using 
Bonfcrroni correction. The analysis of each response variable was carried out using 
the methodology developed in this paper. Based on the above discussion, for each of 
the 6 response variables, we tested the following hypotheses, where i = 1 represents 
COX-l-d, i = 2 represents WT, and i = 3 represents COX-2-d: 



Ho ■ — #2.j — #3j, j ' = lj • • • >4. 

H a :e l!j <e 2>j <e SJ , ., i i. 



Order restrictions on rows and columns of a matrix 



75 



We computed p- values by bootstrapping the test statistic with 50000 resamplings, 
each time sampling the entire record of an animal to preserve the underlying de- 
pendence structure between the 6 response variables. 

After performing Bonferroni correction for 6 tests, significant genotype trends 
were found in subepidermal microblister formation (p = 0.0310), ulceration (p = 
0.0077), epidermal necrosis (p = 0.0034), and acute inflammation (p = 0.0247). We 
failed to see significant trends in hemorrhage (p = 0.1181) and in dermal necrosis 
(p = 0.5170). 

6. Discussion 

In this article we have extended the iterative algorithm of Dykstra and Robertson 
[7] to arbitrary order restrictions on the rows and columns of a matrix as long as 
each row is subject to the same order restriction and each column is subject to the 
same order restriction. However, the order restrictions on the rows need not be same 
as those on the columns. Within each row/column the new algorithm makes use 
of the same estimation procedure introduced in Hwang and Peddada [11]. If rows 
and columns are subject to a simple order then, for independently and normally 
distributed data, the proposed algorithm is identical to the algorithm of Dykstra 
and Robertson [7] . We derive a sufficient condition for the algorithm to converge in 
a single application of Hwang and Peddada [11] method on rows and on columns. 
As an example, the sufficient condition is satisfied in a balanced design. 

Using the point estimators derived in this paper, we introduced a new test statis- 
tic which is a Kolmogorov type distance on the graph of order restrictions as in 
(Peddada et al. [17]). The new procedure is easy to implement and simulations per- 
formed in this paper suggest that it has higher power in most cases than some of 
the existing procedures. A part of the reason for the new procedure to perform bet- 
ter than some of the standard procedures, such as the Kolmogorov-Smirnov test, is 
because it uses improved order-restricted point estimators. As seen from our simula- 
tion studies, the new estimator enjoys substantially smaller quadratic (and quartic) 
risk than the unrestricted estimator, which is used in the Kolmogorov-Smirnov test. 

An added advantage of the proposed methodology is that it is applicable to a 
very broad collection of matrix order restrictions and is computationally easy to 
implement. 
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